Full capacitance-matrix effects in driven Josephson-junction arrays 
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We study the dynamic response to external currents of 
periodic arrays of Josephson junctions, in a resistively ca- 
pacitively shunted junction (RCSJ) model, including full 
capacitance-matrix effects. We define and study three dif- 
ferent models of the capacitance matrix Cp^pr. Model A in- 
cludes only mutual capacitances; Model B includes mutual 
and self capacitances, leading to exponential screening of the 
electrostatic fields; Model C includes a dense matrix Cpyi 
that is constructed approximately from superposition of an 
exact analytic solution for the capacitance between two disks 
of finite radius and thickness. In the latter case the electro- 
static fields decay algebraically. For comparison, we have also 
evaluated the full capacitance matrix using the MIT fastcap 
algorithm, good for small lattices, as well as a corresponding 
continuum effective-medium analytic evaluation of a finite- 
voltage disk inside a zero-potential plane. In all cases the 
effective Cp^pi decays algebraically with distance, with differ- 
ent powers. We have then calculated current-voltage char- 
acteristics for DC-I-AC currents for all models. We find that 
there are novel giant capacitive fractional steps in the I-V's 
for Models B and C, strongly dependent on the amount of 
screening involved. We find that these fractional steps are 
quantized in units inversely proportional to the lattice sizes 
and depend on the properties of C^,^/. We also show that the 
capacitive steps are not related to vortex oscillations but to 
localized screened phase-locking of a few rows in the lattice. 
The possible experimental relevance of these results is also 
discussed. 



I. INTRODUCTION 

There has been considerable recent interest in the 
study of the dynamic response of two-dimensional 
Josephson-junction arrays (JJA) under different phys- 
ical conditions. This interest has both theoretical 
and experimental motivations^. Experimentally, recent 
advances in photoHthographic micro-fabrication tech- 
niques have allowed the manufacture of these arrays 
with specific tailor-made propertieaj. The arrays can 
be made with proximity-effect, i.e., superconducting- 
-metal-superconducting (SNS) junctions, or with 
superconducting-insulating-superconducting (SIS) junc- 
tions. In the SNS case there is essentially zero capaci- 
tance between the superconducting islands. In the SIS 
case the electrodes that form the junctions have nonzero 



self and mutual capacitances that can, as discussed in 
this paper, significantly influence the physical properties 
of the arrays. Theoretically, the arrays are described 
in terms of a large set of coupled, non-linear differen- 
tial equations, overdamped in the SNS case and under- 
damped in the SIS case. The SNS arrays have been stud- 
ied more extensively, for they can be fabricated with more 
ease, better uniformity and larger lattice sizes. In the 
SNS arrays the critical currents can be large and thus 
self- induced magnetic fields must be taken into account, 
through the inclusion of Faraday's-induction law to fully 
describe the experimental systemsoQ. In contrast, in SIS 
arrays the critical currents are relatively low and self- 
induced fields are usually not significant. In this paper 
we consider the SIS case where the electrostatic proper- 
ties of the junctions in the array are of importance, while 
the inductances are not included here. 

When the arrays are driven by DC-I-AC currents, non- 
equilibrium stationary coherent oscillatory vortex states 
may appear, both in the SNS and SIS arrays. Experi- 
mentally, giant Shapiro steps (GSS) have been jpbserved 
in the I-V characteristics of SNS arrays in zeroQ and ra- 
tional magnetic field frustrationsn. The frustration is 
defined by / = $/$o, where $ is the average apphed 
magnetic fiux per plaquette, and $o = h/2e is the mag- 
netic flux quantum. Coherent oscillations of ground-state 
field-induced vortices are responsible for the existence of 
the frcudional GSS when / = p/q^ with p and q relative 
primesQ. This interpretation was successfully verified in 
numerical simulations when the junctions in the array 
were Modeled by the resistively shunted junction (RSJ) 
modefl'B. Half-integer steps have also been found in zero- 
field topologically disordered SNS arrays, due to a special 
oscillatory vortex pattern teemed the axisymmetric co- 
herent vortex state, or ACVSu. Fundamentally different 
half-integer Shapiro steps have also been seen in capac- 
itive array models, with f — 0, when the capacitance 
matrix includes only th* mutual capacitance between su- 
perconducting islandgl3. In all the cases discussed above 
the steps found in the I-V characteristics were microscop- 
ically due to the collective coherent oscillatory motion of 
dynamically or magnetically induced vortices. 

Experimental papers on SIS arrays do report on the im- 
portance of electric field screening, or capacitive effectsi, 
in particular when the junctions are of submicron sizetil, 
but also when they are largeE3. In the cases when the 
electrostatic screening length is smaller than the array 
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size, the screening may affect the array's transport prop- 
erties in a significant way. Most work that has included 
capacitive effects, however, has done so considering only 
nearest-neighbor mutual capacitances. To be more spe- 
cific, consider an array formed by N conducting islands 
where the capacitance of the whole array can be char- 
acterized by a matrix Cp^pi , where f and r' are vectors 
denoting the locations of two electrostatically interacting 
islands. The dynamical equations of motion for a driven 
array that includes general capacitive affects can be writ- 
ten as an extension of the resistively capacitively shunted 
junction (RCSJ) model. This model, which includes mu- 
tual capacitances, has been successful in explaining the 
giant Shapiro steps described above. The corresponding 
equations read: 

13c J2 ^r,r'Or-' + '^p)"^^' + ^ sin {Op+f, - Op) 

r' r' ft 

= ic^t{r.,t) = idc + iacS\n{2TTvt). (1) 

Here Op is the phase of the Ginsburg-Landau order pa- 
rameter for the r-th island; G^^ is the inverse lattice 
Green function in two dimensions (2D) (i.e., the discrete 
Laplacian); the indicates a sum over nearest neigh- 
bors. icxt{r,t) is the external current injected at each 
lattice site r, with both a DC and an AC component that 
oscillates at frequency v. The currents are expressed in 
units of the junction critical current Ic', time is measured 
in units of the characteristic time 1/lUc = h/ (2ei?Jc), with 
R the normal-state junction resistance, e the electronic 
charge, and h = h/2n, with h Planck's constant. The 
capacitance matrix entries are normalized by the single- 
junction dissipation, or Stewart-McCumber, parameter 
(3c — 2eR^IcCm/h, where Cm is the mutual capacitance 
of a single junction. 

The purpose of this paper is to study solutions to 
Eq. (^) for different ever-more-realistic approximations 
to th&.JOf^pi matrix. Just as was done in the inductive 
caseodllj, here we consider the array response to exter- 
nal currents when Cp^p goes from being short- to long- 
ranged. We concentrate mostly on the regime that does 
not show chaotic solutions in the single-junction case, for 
extra complications may arise in that case that can com- 
plicate the analysis further. We also only consider the 
semiclassical regime, where the Josephson energy domi- 
nates the charging energy. As in previous studies we do 
find giant steps in the I-V characteristics in the models 
considered here. In contrast, however, to the GSS de- 
scribed above, the ones described in this paper are not 
associated with collective vortex oscillations. Instead, 
depending on the specific model for the capacitance ma- 
trix, they are related to the phase-locking of specific row 
of junctions to the external drive. 

Here we introduce and study different approximate 
models for the capacitance matrix. We know that we 
cannot exactly calculate the general Cp^p matrix by sim- 
ply giving the geometric configurations of the conduc- 
tors. We start then by assuming that the Cp^pi has only 



the nearest-neighbor mutual component, which reduces 
Cp^p' to a tridiagonal form. This is the model that has 
been studied most in the past (see Ref. Including 
the self-capacitances, which can be done experimentally 
by putting a ground plane underneath the array, leads 
to screening in Cp_pi. The self- plus mutual-capacitance 
matrix approximation reduces the complexity of the full 
Cp^pr significantly. Making this approximation physically 
means that the electric field lines between charges will 
be confined to the plane of the array, with logarithmic 
interaction with distance. Here we want to go beyond 
these approximations and take into account the fact that 
the field lines are experimentally three-dimensional in na- 
ture, so we need to consider further elements in Cp^p' . We 
must then resort to approximate representations of this 
matrix. We have followed several different ways to esti- 
mate the behavior of Cp,pi, some analytic, some numeric. 
On the analytic side, we evaluated the capacitanpe_af 
two disks of finite thickness and radius in a planeoEj. 
We use the result obtained from this calculation to fill 
in, using approximate superposition, what would be the 
matrix elements of a full Cp^pi matrix. This approach 
is similar to the one followed in the fulLinductance cal- 
culations that gave rather good resultsEj. To ascertain 
the nature of this approximation, next we used the MIT 
FASTCAP algorithn£3 that numerically evaluates the Cp^pi 
for small systems. To further understand the fastcap 
results we also evaluated the capacitance of a finite poten- 
tial disk of finite radius embedded in an infinite grounded 
plane. We looked at this problem since it represents a 
type of effective-medium approximation of the fastcap 
algorithm. In all these estimates of the Cp^pi matrix we 
found that it decays algebraically with distance. The rate 
of decay, i.e., the power of the decay is different in the 
various model approximations. This leads to some quan- 
titative difference in the corresponding I-V results, but 
we expect that qualitatively they are correct, and cer- 
tainly when compared with the exponentially decaying 
behavior of the self-mutual approximation to Cp^p. 

In studying the behavior of the I-V's of the full- 
capacitance model, we have found other types of giant 
capacitive voltage steps in the I-V curves, which are of 
a very different origin from the ones studied previously. 
In this paper, we define a row as a series of neighboring 
junctions along the direction transverse to the current. 
Under certain conditions, we have observed these giant 
capacitive fractional steps (GCFS), which are caused by 
phase-locking of junctions in a given row to the exter- 
nal drive, with no significant phase-locking between rows. 
The array exhibits the phase-locking behavior which un- 
derlies the GCFS, but we find that this does not occur 
throughout the entire array. Rather, it only occurs in 
certain rows near the edge of the array, due to screening. 
We discuss these results in more detail in the main body 
of the paper. 

The outline of the paper is the following: In Section II 
we present the equations that describe the different mod- 
els mentioned above, as well as briefly discussing the 
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methods we used to study them. We also present some 
analytical results which we used in our analysis and to fill 
out the entries of the full capacitance matrix. Section III 
discusses the bulk of our results and their analysis, for 
the cases where the driving current is purely DC or a 
combination DC+AC. We study the I-V's, spectral func- 
tions, and topological phase and current distributions, to 
support our view that the steps are indeed generated by 
phase-locking of a few rows of junctions (rather than the 
entire array), and not by the presence of a stable vortex 
configuration, as is the case with other giant fractional 
steps. Finally, in Section IV we summarize and discuss 
our results. 



II. CAPACITANCE-MATRIX MODELS 

In this section we proceed to define the different capac- 
itance matrix models studied and compared in this pa- 
per. Our goal is to solve Eq. (|l]) for different Cp^p matrix 
models. We begin by converting the set oi N = x Ly 
second-order ordinary differential equations into 2N first- 
order equations. Explicitly, we define the variables zp and 
zp that satisfy the equations. 
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where the indices r, r', and r" run from 1 to N, that 
is, over the whole array. Periodic boundary conditions 
are used in the direction perpendicular to the current. A 
schematic of the arrays under study is shown in Fig. |l|. In 
order to simulate large arrays, it is useful to use the gen- 
eral properties of the capacitance matrix that arise from 
the positivity of the total electrostatic energy of the array 
of electrodes. Specifically, Cp^p' > 0,V r = f , Cp^p' < 0, 
\/ f ^ r' , and Cp^p > J2r'^f\^f,f'\- In addition to these 
properties, we have translational invariance in a periodic 
array, so that we can write C{r,r') = C(\f— r'j), which 
allows for a significant simplification in the calculation 
of Cp^pi. Our approach parallels the scheme used in the 
dynamic study of inductive Josephson-j unction array^, 
and we consider here the following three distinct capaci- 
tive models: 

Model A: In this case we consider only the contribu- 
tions by the mutual capacitances between electrodes in 
the array. Here the capacitance matrix is simply tridiag- 
onal, i.e., 



r,r'±e^ 



r,r' ±ey 



^-1 



with Bx and Cy are unitary vectors along the x- and 
y-directions. This model has been studied by many 



authorffflli, both for ordered and disordered arrays, 
and some results are already known. 

Model B: This model extends Model A by including the 
self-capacitance term in the capacitance matrix as: 



where rc = \Cs/Cm\^ and Cs and Cm are the self- 
capacitance and the mutual capacitances, respectively. 
This model tends to Model A in the limit Tc 0. In the 
continuum limit the inverse of the capacitance matrix 
becomes 
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with Kq{z) is the modified Bessel function of zeroth or- 
der. Asymptotically, we have that Ko(z) ^ —£nz, for 
z <^ 1, and Ko{z) for z » 1. We see that 

when the self-capacitance is very small the interactions 
are essentially logarithmic with distance, and for large 
Cs they are exponential. We define then the screening 
length as Aq — We will describe our results for the 
models studied here as functions of the parameters, /3c, 
v and Aq, or Vc- 

Model C: Here we consider the capacitive effects of all 
conductors on each other. This means that the capac- 
itance matrix is dense, and the specific form of the in- 
dividual entries depends on the details of how we model 
the conductors in the array. We have used several routes 
to analyze the long-range interaction that gives approx- 
imate full capacitance-matrix models. Most of our I-V 
calculations are based on a Cp^pi obtained by an approxi- 
mate superposition model from an exact analytic solution 
for the capacitance of two circular disks of finite thickness 
and radiusoE3. We review the derivation of this result 
in detail in the Appendix. 

We place each disk at a lattice site in a square array. 
For the total capacitance matrix we use simple superpo- 
sition of the capacitances between any two disks, as done 
in the inductive case, so as to get 



Cp^p' — C{i, j) 



4 tanh 



1 - 



and 



C'rr 



C(i = l,j-0) 



4 tanh" ^ (VI -D^) 



(4) 



(5) 



where f — f = iBx + jfiy, i and j are integers, D is the 
diameter of a disk and R is the distance between their 
centers (see Fig. ||). Note that the magnitude of D defines 
the packing of the disks as shown in the figure. The two- 
disk expression for C{i,j) given above allows us to go 
from weak to strong screening by just varying D. For 
different values of D, C{i,j) has an initial rapid decay 
with distance and then it decays algebraically. We show 
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in Fig. H the behavior of the effective full C{i,j) matrix 
for some values of D. 

For comparison to the two-disk model result, we also 
have used the FASTCAR^umerical capacitance extraction 
tool developed at Mlllla. This algorithm quantitatively 
calculates the capacitance matrix for a given distribution 
of conductors. Whan et al. [ |l^ have studied one and 
two-dimensional arrays of circular dots using fastcap. 
They also obtained a power-law decay with exponent 
close to one. We also use their result in the evaluation of 
DC IV's in Section III. We have also explicitly used FAST- 
CAP, to compute Cf^pi for an 8 x 8 array of square plates 
of various thicknesses, ranging from 1000 A to 35000 A. 
We took the dimensions used by van der Zant et al&3. 
A log-log plot of the capacitance versus the distance is 
shown in Fig. ^. This capacitance result tends to the 
two-disk array model for a specific value of D, but differs 
for others. The essence of the result, however, is that the 
off-diagonal capacitance matrix decays algebraically. 

The FASTCAP approach is limited in that the arrays 
that can be simulated are relatively small. (For example, 
11 X 11 was the largest size possible on a fast DEC Al- 
pha with 196 MB RAM.) In the FASTCAP algorithm the 
capacitance matrix entries are calculated iteratively, by 
successively grounding all but one of the conductors, and 
holding it at a fixed nonzero potential. The limitations 
of this approach are that it is time-consuming and not 
many parameter values can be easily considered. In con- 
trast, the two-disk model is more flexible and allows the 
study of larger lattices and a wider range of parameter 
values. 

In order to further understand how the capacitance 
matrix behaves we have also analytically calculated the 
electrostatic potential, and then the capacitance, of a cir- 
cular disk of finite width set at a fixed potential embed- 
ded in a conducting plane set at zero potential The disk is 
of radius s = D/2 located at the coordinate origin and at 
a finite nonzero potential Vq , while the conducting plane 
is at zero potential. We considered this problem since 
it is, in a sense, similar to the lattice algorithm used in 
FASTCAP and it also gives us another evaluation of the 
capacitance in the continuum limit. We want to calculate 
the potential outside of the plane and thus start from the 
general expression for the potential at point r given by 
(see Jacksoro) 



dA'. 



Here F{r,r') is the Green function for two unit charges 
at z and 2', n is the normal vector to the surface of inte- 
gration S, and dA' the surface differential. In cylindrical 
coordinates (p, ip, z), F{r, f") is given by 



F(p,p',0,0',z,z') 



1 



\/ + p'^ — 'i.pp' cos 7 + (z — 
1 

p^ + p'-^ — 2pp' cos 7 -I- (z -I- z')2 ' 



where 7 = [ip — ip'). After evaluation of the z-derivative 
of F we get 



zVop'dp'dip' 



(p2 _|_ pa _ 2pp' cos 7 + z2)3/2 ■ 
The surface charge is given by the expression 
1 dV{p,ip,z) 



47r 

which gives the integral. 



dz 



z=0 



p' dp' dip' 



/o JO [p'+p"'-2pp'co8ip-p')]y^- 
The yo'-integral can be evaluated explicitly giving 




sin^{p — + — 2ps cos{p — p') 



27r 



cos((p — p')dp' 



P Ja sm'^{p — p)')\J p^ + — 2pscos{ip — p') I 

The first integral is zero while the last two integrals are 
elliptic in nature and are given in tables. Here we are 
interested in the long-distance behavior and thus expand 
the integrals in the index of the elliptic integrals. We 
can then obtain the corresponding capacitance, defined 
as charge divided by voltage and the volume of the disk 
of thickness h, as 
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Sn^h \ p^s 



6/s2(4p2 _ 2s2 _ 2sp) 
ip-sfip^-s^) 

8p6s4(6p2„4g2„25p) 
(p-s)13(p2-s2) 



(p_,)4(^2„,2) 



(6) 



One of the advantages of this result is that it shows that 
there are different algebraic contributions depending on 
the p-range. We show in Fig. ^ that for small p this 
expression is almost linear in the log-log plot but not 
exactly so for larger distances. 

All the results for the full capacitance analyses are 
shown in Fig. ^. There we see that the rates of decay 
are different depending on the approximation but the ba- 
sic nature of the result is the algebraic or nonexponential 
decay of the off-diagonal elements of the capacitance ma- 
trix. Note, however, that the comparisons between ca- 
pacitive models here are qualitative since the definition 
of the parameters in the different models is not exactly 
the same. We conclude that the two-disk model contains 
the essence of a full capacitance matrix and thus have 
done most of our calculations using this model. 



4 



A. Calculational approach 

In solving the dynamical equations of motion a good 
amount of time is spent in calculating the inverse of the 
capacitance matrix. Inversion of Cp^pi may be tackled in 
different ways, depending on how much we know about 
it. 

In the case when Cp^p is tri-diagonal, it has only mu- 
tual or mutual+self capacitances and its Fourier repre- 
sentation in momentum space is also tri-diagonal. Taking 
advantage of this form, we can use an^ptimized (lower- 
upper) LU decomposition algorithnJij to invert Cr.r' , 
such that the calculation time grows like ~ iVlogA^, 
rather than ~ N'^ . In the case when Cp^pi is dense, we can 
still use the translational invariance to go to momentum 
space, although the multiplication (again by LU decom- 
position) will now grow like ~ iV^. 

Another possibility is when we already know the in- 
verse of Cp^pi and we just need an efficient way to carry 
out the multiplications. Once again, translational in- 
variance can be used to reduce the complexity of this 
task. We can rewrite equation Eq. (^) using the defini- 
tion s{r) = w(r) - Y.r' Gp^pZp - sin (6'^^+/i - Bp), 
as 



1 

Jc 



(7) 



As pointed out by Phillips et aln, this has the form of a 
convolution, and can thus be evaluated in ^ N"^ multipli- 
cations, plus the cost of performing the Fourier transform 
(~ TV log N) , as opposed to ^ for the direct multipli- 
cation. Because we do not have periodicity in the-jcurrent 
direction, in order to obtain a linear convolutionl3 we use 
the standard 'zero-padding' techniques as is done in dig- 
ital signal processing. The integration in time is carried 
out mostly by the second-order Runge-Kutta (RK) algo- 
rithm. 
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and the current vorticity by 



c(i?) = J2 sin (A6I) 

P{R} 



(9) 



(10) 



Here P{R} denotes the sum around the four bonds en- 
closing the plaquette at position R, and nint(x) is the 
nearest integer to x. It is important to remember that 
the phase difference A9 is restricted to the range [— tt, tt). 
The fractional steps which we have found in the capaci- 
tive models are not vortex related but it was important 
to calculate the vorticities to check that. We do not show 
figures since there is nothing to see. 

In particular in Model C we find that there are neither 
topological vortices nor eddy current vortex produced in 
the fractional steps. This appears to be true even when 
we let the system evolve for 500 periods of the driving 
current, long after the voltage has stabilized. In addi- 
tion, in the FASTCAP model, when we look at the eddy 
vortices, we find that, apart from a short-lived transient, 
no vortices are present. 

To further understand what is happening we also have 
looked at the spectral function defined as 



S{h') — lim 



(11) 



For integer steps S{iy) has peaks at frequencies nv, for 
periodic and disordered arraysQ. This is consistent with 
the resonant nature of the integer steps. For the dis- 
ordered arrays that produce the ACVS state there are 
half-integer vortex steps, corresponding to an S'(z^) that 
has additional peaks at frequencies nv/2. 



B. Physical quantities calculated 

One of the quantities we are interested in computing 
is the total voltage drop per junction in the array given 
by the Josephson relation, 

<""»^£mi^ I (""■-' -''")^ 

The Josephson-j unction arrays studied before exhibited 
steps of two very different origins: the integer steps, due 
to coherent phase-slippage, and the fractional steps, due 
to vortex motion. It is then important to see whether 
or not there are vortices associated with possible quan- 
tized fractional steps in the I-V's. We define two types of 
dual lattice vorticities: one topological and one current 
related. The topological vorticity is defined as 



III. RESULTS 

In this Section we consider the response of Models A-C 
under different current-driven conditions. First we con- 
sider the case when the array is subjected to a DC ex- 
ternal current and next when there is a DC-I-AC current 
drive. We also look at other physical measures defined 
above, that further help us in the analysis of the results. 

A. DC current— voltage characteristics 

Josephson-junction arrays driven only by a DC current 
have been studied before in the overdamped case and also 
in Model vlld. Here we examine the I-V's for 8 x 8 arrays 
in Model C based on the two-disk approximation as well 
as a square plate model based on a FASTCAP calculation. 
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We have computed the DC I-V's for Model C in a DC 
drive for different values of the parameter, D. When 
D < 1, and with Aq ^ Lx,Ly, there is basically no 
scrcoaing and thus the I-V's have the usual quadratic 
formll3: 




In Fig. Q we show the behavior of the I-V's for two ex- 
treme values of the disk diameter D. As D slowly de- 
creases from 1.0, we find little change in the shape of the 
I-V's until D ~ 0.7. At this point, the I-V rather loses 
its smooth curved shape, and assumes a piecewise linear 
form characterized by a very sudden increase, and then 
to linear ohmic behavior as the current is increased. We 
also note a slightly lower critical current for the latter 
case. 

Next we also calculated the I-V in the fastcap disk 
model considered in Ref. ^ in a regime close to that of 
Model C considered above. In this case the results shown 
in Fig. ^ are somewhat more tentative, for we may be 
operating outside the regime in which the fastcap model 
is valid. Qualitatively, the I-V's display the transition 
from a smooth shape for a large value of D (defined as 
a normalized diameter of the disk) to piecewise linear at 
smaller l)-values. Note that quantitatively this length D 
is not exactly the same as the one defined for Model C. 
The conclusion from these calculations is, however, that 
the results evaluated with both models are qualitatively, 
if not quantitatively similar. 

B. DC+AC current— voltage characteristics 

Here we present results for DC-I-AC driven capacitive 
JJA and compare results for Models A, B, and C. Model 
A has been studied to some extent and we briefly mention 
those results here, but we will be mostly concerned with 
results for Models B and C. 



1. Model A: Mutual capacitances only. 

This case has been studied in Ref. |l^ mostly in the 
transition region between regular and chaotic behavior 
for the single-junction problem. Ajiart from the integer 
giant Shapiro steps Hagenaars et alS3 find the formation 
of a half- integer step, within a range of values of Pc which 
include chaos in the single junction, with characteristics 
similar to those of the ACVS state mentioned in the in- 
troduction. In our analysis we mostly stay away from 
the single-junction chaotic state, and do not see the half- 
integer steps. Next we look at what happens when there 
is screening in the capacitive models. 



2. Model B: Self+mutual capacitances only 

Here we present results when the capacitance ma- 
trix is composed of the self capacitance of each island 
to ground plus the mutual capacitance of each island 
to only its nearest neighbors. The time-step used in 
the second-order Runge-Kutta calculations was 0.05 for 
all the results in this section, with a current grid of 
Sidc/ic — 0.005. Except where otherwise noted, averag- 
ing was carried out over 6x 10^ time-steps, after throwing 
out the first 2 x 10^ for equilibration. 

As mentioned before, addition of the self capacitance 
to Model A introduces exponential screening that for 
large distances is measured by the screening length Aq — 
\/ Cm/Cs — 1/y/rZ- Our results were obtained princi- 
pally on an 8 X 8 array, but we have seen similar behav- 
ior in much larger arrays (e.g., 40 x 40, see Fig. ^). We 
find that there are new giant capacitive fractional steps 
(GCFS) as we vary Aq. When such steps occur, they ap- 
pear as multiples of 77^37, where n = 1, 2, 3 and possibly 
4. 

The results obtained when Aq 3> Ly were, as must be, 
identical to those obtained with the capacitance matrix 
composed of Model A. The existence of these integer- 
steps depends on the locked-in coherent phase oscillation 
of all junctions in the array with the drive. On reduc- 
ing the screening length to a value smaller than the ar- 
ray size, the integer-steps began to disappear, and then 
new steps appeared at fractional values of the normalized 
voltage. 

It can be seen in Fig. ^ that, for Aq > Ly, the width 
of the integer-step is entirely unaffected by screening. As 
Aq approaches ij,/2, the step at n = 1 begins to shrink, 
and as Aq ^ 1, the step is rapidly destroyed. In Fig. ^ 
we show the change in the width of the n — 1 step as a 
function of Aq. 

We then examined the step width for the GCFS |-step 
as functions of the frequency (3c and Aq. As can be seen 
from Fig. ^ the appearance of this step is quite strongly 
dependent on the drive frequency, and it is absent for 
frequencies 15% greater than or less than ly = 0.1, with 
/3c = 0.50. Next we fixed v, and looked at the width of 
the fractional steps, as a function of (3c- Our results are 
shown in Fig. |lO[ The j- and |-steps are seen to ap- 
pear for Pc ^ 0.3, reach their maximum values at around 
(3c ~ 0.7, and then slowly disappear, as /3c is increased. 
All fractional and integer steps have disappeared by the 
time (3c has reached the value 3.0. Each step has a max- 
imal screening length denoted A^^^ for which its width 
is maximized, and the lower the order of the step, the 
smaller its A™^'^. We mostly stayed outside the range 
of (3c for which the single-junction behavior is chaotic, 
though Fig. |l^ shows that chaos is not necessary to gen- 
erate the GCFS. 

We then fixed both i/ = 0.1 and (3c — 0.6 (parame- 
ters for which the step- widths appear to be maximized), 
and looked at how the step-width varies as a function of 
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Aq. The results are shown in Fig. |Tl|. We decreased the 
screening length from a large value until Aq had reached 
half the array size, i.e., Aq ^ 4. At this point, the frac- 
tional step size jumps up rather rapidly to almost its 
maximum value, remains there until the screening length 
has reached the value 1.0, and then falls rapidly to zero 
once again. This behavior can be seen to happen for both 
the J- and |-steps. 

We obtained these steps by starting up from a lower 
value of idc, since one cannot see the step by simply set- 
ting the value of idc according to the expected step values 
in the I-V's. The fractional steps are hysteretic, as indi- 
cated by the inset in Fig. |ll|. This is a typical property 
of Model A arrays, and has been seen in previous stud- 
ies, whose authors interpret them as evidence of vortex 
inertiall3. Both of these facts indicate that the GCFS are 
metastable, since their existence depends on the previous 
history conditions of the array. 

Usually, the appearance of a step at a voltage V — - 
indicates that a process with frequency underlies the 
production of the step. This should be visible in the 
Fourier power spectrum of the array voltage. For the 
integer steps that is what we see but not so for the GCFS, 
which signals that they are of a different nature to the 
steps studied before, as we further discuss below. 

3. Model C: Full capacitance matrix 

In this subsection, we present our results for the ca- 
pacitance matrix composed of all nonzero elements. In 
this case, the RK time-step calculation used was 0.05 
for all the results in this section, with a current grid of 
5idc/ic = 0.005. Except where otherwise noted, averag- 
ing was carried out over 6 x 10^ time-steps and the first 
2 X 10** were used for equilibration. Longer time-series 
have shown that while there are some long-lived tran- 
sients, eventually the state becomes stable. 

Since the off-diagonal matrix elements of Cp^pi decay 
algebraically, the definition of screening length as done 
in Model B does not strictly apply. Nonetheless, because 
of the specific decay in Model C we can still quantify the 
screening by the ratio of the off-diagonal capacitances to 
the diagonal ones, using the definition of Model B. We 
could have used another definition in which we sum over 
all the matrix elements, normalized appropriately, but 
the end result, that is the measure of screening, would be 
about the same. (In fact, we did just that, and found that 
while it resulted in a rescaling of the screening length by 
a small numerical factor, it did not change the qualitative 
behavior.) 

Our results were obtained principally on an 8 x 8 ar- 
ray, but we have seen similar behavior in much larger 
arrays (e.g., 40 x 40). In this model, in contrast to Model 
B, we have seen only the |-step which has to do with 
the nonzero decay of the off-diagonal matrix elements of 
the capacitance matrix. In Model C we can move con- 



tinuously from Model A to Model B by adjusting two 
parameters: the self capacitance (by varying Tc) and the 
diameter of the disk, D, expressed in units of the lattice 
spacing. D takes values in the range [0, 1), the upper 
limit giving Model B (or Model A if, in addition, Tc = 0). 
The lower limit corresponds to an array of superconduct- 
ing dots which are too far separated for tunneling. As 
we decrease D from 1, for Aq ^ La, a voltage step ap- 
pears at {V) = |, for _D = 0.997, and disappears for 
D = 0.95. The results are shown in Fig. for Aq = 32 
on an 8 X 8 array. The maximum width of the S-step is 
of the same order of the half-step seen beforcEEl, though 
we chose Pc — 0.5 specifically toie outside the chaotic 
regime for which Hagenaars et aZJl3 saw that step. 

The magnitude of the step appears to be 77-^ , based 
on runs for Ly = 4, 8, and 16: that is, for each lattice size, 
we observe the appearance and subsequent disappearance 
of a step at (V) = 77-^ as D is varied down from 1. 
The occurrence of Ly in the expressions suggests that the 
GCFS are edge-like effects. They do not occur for Model 
j4, but do occur for certain parameter values in Model 
B, as can be seen in Fig. |^. We also found qualitatively 
similar behavior (i.e., the existence of a step only for 
a limited range of Z? ~ 1) using the expressioupJor the 
inverse capacitance matrix model of Whan et alJiB 

We have examined the width of the |-step as a function 
of Aq, as we fixed both ly — 0.1 and f3c = 0.5. The results 
are shown in Fig. |l^. There are some clear differences 
between the behaviors of Models B and C. We stress that 
because the screening in Model B is exponential for long 
and logarithmic for short distances, as compared to the 
algebraic decay in Model C, this is clearly seen in the 
different results in this figure. 

The GCFS steps appear to be related to a partially 
screened phase-locked state. They do not appear to be 
metastable, in the sense that there is no need to 'step' up 
onto them, as was found in Model B. We also produced 
animations that show no stable vortex patterns, but in- 
dicate that all junctions in a transverse row evolve in 
phase without an apparent simple relationship between 
adjacent rows. We examined the voltage states of all 
junctions along the central column of the array (that is, 
the ladder consisting of all the junctions with 2: = 4), and 
found that there are essentially two kinds of behavior: 
The behavior of the edge junctions seems to be symmet- 
rical about the line y = Ly/2, with their voltage (normal- 
ized for frequency) equal to unity. The other junctions, 
while not locked in or very much out of phase with each 
other, are in zero-voltage states. This is clearly due to 
screening. 

We also looked at the spectral function of a Model C 
array. We find, sitting on the |-step (having stepped up 
onto it), that when we study the spectra of the junctions 
within the center row, at one of the edges (all within the 
ladder defined by a; = 4) or the whole array, all three 
cases show similar behavior. 

There is a small but very long-lived transient, which 
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shows up in the power spectral density, as two peaks 
at frequencies approximately equal to ~ ^ and -y- (as 
shown in the inset of Fig. |lj) . Also present in the spec- 
trum are the expected drive frequency v = 0.1, and its 
harmonics. We have verified that it is a transient and is 
not responsible for the appearance of the fractional steps, 
by allowing the array to equilibrate for very long periods 
of time. Under these conditions, we still see the step, and 
the power spectrum is very clean and shows no trace of 
the transient, as shown in the inset of Fig. 

Figure |l^ shows the I-V's of an edge row of junctions, 
an interior row and the average of the whole array. Both 
edge rows behave similarly to each other, but differenatly 
from the interior rows. It is clear that the edge row enters 
a nonzero voltage state at a lower current than the inte- 
rior row (i^c = 0.81 versus idc = 0.885). That we have 
two out of seven rows in a voltage-producing state, and 
five out of seven in a zero- voltage state, seems to indicate 
that the |-step has its origin not in a coherent oscillation 
of the array as a whole, but rather in the fact that some 
of the junctions get driven into a nonzero voltage state at 
a lower current than others, producing the observed |- 
step. This is more easily seen in the 16 x 16 array, where 
multiple steps are visible. In particular, for D = 0.999, 
we see steps at ^ and j^. Looking at the I-V's for 
the individual rows, we can see first the outermost (edge) 
two rows (one at each boundary) switching on, followed 
by the next outermost rows, followed by the next out- 
ermost ones. The points where these rows sequentially 
switch on correspond to the beginning of a new step. 

4- Size dependence 

In the results reported above the size of the system ap- 
pears prominently. The question is then how stable are 
these results as a function of lattice size and parameter 
values. Most of the results described in the previous sec- 
tions have been generated in lattices of size ~ Ly ^ 8. 
The steps we find are multiples of ^ (since Ly — 1 = 7). 
We also have looked at array sizes of = Ly = 16,40. 
In both cases, we found analogous GCFS at integer mul- 
tiples of -j^ and ^ . In Model B, we notice that the maxi- 
mum step width for the -^-step in the 8x8 array is about 
Aide = 0.10, whereas for the 40 x 40 array, we find it to 
be ~ 0.085. The widths are therefore comparable, and 
seem largely independent of lattice size. Figure |^ shows 
the I-V for a 40 x 40 array, with a screening length much 
smaller than the array (Aq = 1.75). 

In Model B we have examined the critical screening 
length A"'* for which the first step appears as a func- 
tion of size Lx- We show in Fig. |l^ the corresponding 
data, along with an approximate fit that gives A"'* ~ 
^o.9i4±o.058^ Figure |l^ shows the maximum width of the 
first step as a function of for a given maximum screen- 
ing length A™^'^. Note that the first step is generally the 
biggest, regardless of lattice size. In addition, we show 



in Fig. |/7| that the maximal width of a step of a given 
order appears to saturate once the lattice size is above a 
certain minimum. 

For Model C, we also investigated the size of the frac- 
tional steps as a function of the disk diameter, D. Fig- 
ure |l8| shows the dependence of the width of the j^- 
step in the 16 x 16 array. The qualitative behavior is 
very similar to that of the 8x8 array. Quantitatively, 
the difference lies in the values at which the steps ap- 
pear and subsequently disappear. In the 8x8 array, 
the region of I?-values for which the steps are visible is 
D £ [0.950,0.995], while in the 16 x 16 array, it is much 
smaller than that: D e [0.992,0.9999]. 

IV. SUMMARY AND DISCUSSION 

In this paper we have considered the general response 
of a Josephson-junction array that incorporates the ca- 
pacitance matrix at different levels of approximation. We 
have considered three different models. Model A includes 
only the mutual capacitance between nearest-neighbor 
islands-. This model has been studied to some extent 
beforeEHl. Model B includes the diagonal terms of Cp^p 
and had not been studied before. The important el- 
ement in Model B is the exponential screening present 
at long distances that introduces important changes in 
the current response in the array: to wit, the new gi- 
ant capacitive fractional steps (GCFS) that arise in the 
I-V's in this model, appearing and disappearing as the 
screening length Aq changes. Model C incorporates the 
electrostatic interaction between all conducting islands. 
Since there is no exact solution known for an exact full 
capacitive matrix, we considered different model approx- 
imations to the full Cf^pi. Most of our results were based 
on an explicit solution for the capacitance of two disks 
of finite size in the plane using a superposition to con- 
struct the full matrix. This approach is similar to the cme 
previously used in the full inductance-matrix problemEj. 
We also used the |iuimerically exact fastcap algorithm 
developed at MITll3. The essence and main result from 
all these calculations is that the full capacitance decays 
algebraically rather than exponentially as in Model B. 
We then concluded that the two-disk approximation con- 
tains, qualitatively at least, the essential properties of the 
full capacitance matrix of the more accurate problem. 

Our results indicate that for screening lengths Aq much 
larger than the lattice size, the TV's display the char- 
acteristic giant Shappji steps of the mutual capacitance 
model only {Model j4E2l). That is, there are giant integer 
steps, and also (for certain values of f3c) the much smaller 
half-steps that are triggered by the onset of chaos in a 
single junction. The half-steps display the characteristic 
ACVS. The essential new result from our calculations in 
Models B and C is that there are new GCFS that essen- 
tially depend on the amount of screening in the model. 
The GCFS are not produced by oscillating vortices, as 
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in the non-fuU-capacitive problems, but are localized in 
restricted areas of the lattice due to electrostatic screen- 
ing. 

Generally, the size of the GCFS varies inversely with 
(though not necessarily proportionally to) the order: that 
is, lower-order steps are larger than higher-order ones. 
These fractional steps in Model B are metastable in na- 
ture, (that is, history-dependent but stable once pro- 
duced) but remain for the longest runs we have carried 
out. Once Ag ~ 1, all steps begin to disappear, and by 
Aq = 0.5 they have completely vanished. 

For Model C when D = 1.0 (and for large enough 
screening length, Aq), we have the under-damped results 
of Model A, and only see integer steps in the I-V's. As D 
decreases from 1.0, we enter a regime in which different 
rows have different values of i^c at which they assume the 
coherent, nonzero voltage state. The rows closer to the 
boundaries of the array maintain the critical current of 
the under-damped array, while those towards the center 
become increasingly difficult to switch on. The inter- 
play between the switching-on values for different rows 
determines the relative sizes of GCFS that we see. There 
may be no value of idc for which all the junctions in the 
array can oscillate in phase. This corresponds to both 
the shrinking of the integer steps, and the appearance of 
the fractional steps. As D decreases further, we enter a 
regime in which the coherent oscillation is no longer pos- 
sible even for a single row. At this point, we have lost all 
the steps, fractional and otherwise. 

The question naturally arises: To what extent does 
this behavior depend on the model used to construct the 
capacitance matrix? We believe this behavior to be in- 
dependent of the model used, based on simulations we 
have carried out using results obtained independently, 
and by a completely different method by Whan et alx^ 
and the fastcap algorithm calculations. The idea now 
is to check the predictions of this paper experimentally. 
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APPENDIX: A 

In this appendix, we describe how to obtain the ca- 
pacitance of two disks given in Eqs. (0) and (^) used in 
the array calculations. It is best to use bipolar coordi- 
nates (see Ref. ^0|) to obtain the relevant expression but 
then we shall transfppn to Cartesian coordinates as used 
in the calculationsEjH. These coordinates are related to 



a sinh 77 



a sin ^ 



the Cartesian coordinates through the equations 



(cosh 77 — cos^) ' 



(cosh?7 — cos^) ' 



(Al) 



The parameter rj measures the diameter of the disks. 77 = 
—00 represents an infinitesimally small disk at position 
(— a,0). As ?7 increases, the disk grows becoming infinite 
in size, at 77 = 0. Increasing 77 further causes the disk to 
shrink until it becomes a point at (a, 0), for 77 — +00. The 
parameter a denotes the distance of the points rj = ±00 
from the y-axis, a coth 77 is the distance from the center of 
the disks to the y-axis, i.e., half of the separation distance 
between the centers of the disks. The lines of constant 
represent the boundaries of the disks. 

The orthogonality of the bipolar coordinates can be 
establishe d by evaluating the infinitesimal displacement 
from Eq. (Al) and identifying the respective scale factors 
and unit vectors: 

dr = ixdx + iydy + Czdz = e^h^d^ + erjhridrj + Czhzdz, 



h. 



1, 



cosh 77 — cos ^ ' 
- Cx sinh T] sin ^ + e j, (cosh 77 cos ^ — 1 ) 

cosh T] — cos f 
- e^ (cosh T] cos ^ ~ I) — Cy sinh 77 sin ^ 



cosh 77 — cos ^ 

The electrostatic potential function between the two 
disks satisfies the Laplace equation, so from the previ- 
ous equations we get, 



1 



h^hrj 



</.(^, 77)^0. 



(A2) 



The general solution to Eq. (A2) can be written as 



</'(C: V) = ^ [^m cos(77i^) -|- B„i sin(r7T,f )] 



rn—0 



X [Cm cosh{mi]) + Dm sinh(77277)] . (A3) 

Taking the disks with centers at 77 = 771 and r; = 
772, the boundary conditions are (/>(^,7y = 771) = Vi, and 
(/)(^, 77 = 772) = ^2 = 0. To get the electrostatic potential 
function we only need the case m = 0: 



V2{r] - 771) + Vi(?72 - 77) 



m - m 



(A4) 



The electric field intensity is obtained by taking the nega- 
tive gradient of Eq. ( A4), and only has e,,-direction com- 
ponents: 



V2-V1 



(A5) 



The electric charge distribution on each of the disks is 
given by 
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Brj ■ E 



An 

{V2 - V^i) (cosh 771 - cos^) 



ATr{T]2 - rii)a 



in ■ E 



V=V2 

(V2 - Vi) (cosh 772 - cos^) 
47r(?72 - rii)a 



(A6) 



The total charges are obtained by integrating Eqs. ( A6) 
over the respective disks 



h i'2-iT 

/I = / / cr(C,7y = iji)h^d^dz 
10 Jo 

V2~Vi 



47r(772 - 771) 



-27r/i : 



V2-V1 



2(?72 - m) 



(A7) 



Therefore the capacitance of the bipolar capacitor with 
height /i, follows from Eq. (A7), 



Qi 



{Vi - V2)h 2{t^2 - m) ■ 



(A8) 



Now to use this result in our calculations we need to 
transform it appropriately to the lattice structure of the 
square lattices considered in this paper. Eq. (AS) as- 
sumes that the two disks are different. In the arrays we 
have studied the disks have the same radius and thus 
r]i ~ 1] and 772 = ^Vi = so that the capacitance of 
the two disks is 



C 



1 

Ar]' 



(A9) 



Now going back to Cartesian coordinates, we can express 
ry and ^ as follows: 



77 = tanh 



^ = tanh 



2ax 



2ax 



? ? 9 

_ _ yZ 



(AlO) 
(All) 



We define d as the ratio of the radius of the disk to the 
distance of its center from the origin (or, equivalently, the 
ratio of the diameter to the distance between the centers 
of two adjacent disks): 



d = 



1 



a|csch77| 
a coth 77 cosh 77 



(A12) 



This provides a way to parameterize the separation be- 
tween adjacent disks. Since coshx = (e"^ + e^^)/2 > 
1, Vx, d lies in the range [0, 1). When d < 1, the disks 
are almost touching, and each disk occupies almost all 
of the lattice unit cell. We expect this to correspond to 
Model A. When d ~ 0, then the disks are separated by 
a very large distance in relation to their diameter. This 



scenario corresponds to weak screening. We will express 
the capacitance in terms of this d parameter. 

Taking y — 0, an d measuring distances in units of a 
hereafter, Eq. (AlC) becomes 



1] =^ tanh ^ 



2x 



1 



(A13) 



The points where the disk intersects with the x-axis are 
given by {x — coth 77)^ — csch^77 (the equation of a circle 
of radius csch?;, centered at coth 77). This has solutions 
x = (cosh 7/ ± l)/sinh77, which can be expressed in terms 
of d: 



l±d 



(A14) 



(There are two values of x for each d — each gives the same 
expression for the capacitance, as it should.) Plugging 
the solutions in Eq. (A14) into Eq. (A13) gives us: 



rj = tanh ^ 

= tanh^"'^ 



'2(1 ±d) 
\/\~d^ . 



l-rf2 



\\±df 



Therefore, the capacitance behaves as: 



4tanh~i [VI - d^\ ' 



(A15) 



(A16) 



where d is the ratio of the disk diameter to the distance 
between the centers of two adjacent disks. This is the 
result used in our calculations and given in Eqs. (^ and 

(!)• 

As a check on this result, we can look at the behavipc 
as d ^ 0. This problem has been solved by Jacksonta 
(problem 1.7, p. 51) who gives the capacitance per unit 
length for a pair of infinite parallel wires, with a ratio of 
radius to separation d <C[1, as C = —1/4.^7^^. 

By using the expressiontll: 



tanh ^ z = —In 



1 



1 



and expanding binomially the square root of Eq. (A15) 
in the limit of small d, we find that 



tanh 



\4n = in {2/d) ~ -In d. 

(A17) 



This gives us 



C{d) 



-1 



Atnd 



(A18) 



Thus, we have agreement between our model and a 
known solution in the limit of small d. 
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FIG. 1. Diagram of the array studied. The symbols Kl de- 
note the capacitive junctions. Current is injected from the 

bottom, and extracted from the top, with periodic boundary 
conditions along the direction perpendicular to the current 
flow. 

FIG. 2. This figure illustrates the two-disk model. D is 
their diameter measured in units of the lattice spacing and 
R the separation distance between them. When D < 1, the 
array is close-packed, with D almost equal to the lattice spac- 
ing. When D <^ 1, the array is loosely packed, with the disks 
widely spaced from each other, and thus no current can flow 
between them. Our model is not valid for this value of -D ^ 1. 

FIG. 3. Comparison of behavior of capacitances for dif- 
ferent model calculations. The lines represent the following 
models: FASTCAP calculation ( ); Model C two-disk ap- 
proximation ( ) for D = 0.99, and ( ) for D = 0.1; 

flnite-potential disk embedded in an infinite zero-potential 
conducting plane ( ). The inset shows the conduc- 
tor plate array used in the FASTCAP calculations. The plates 
are of side £, with thicknesses ranging from 1000 to 35000 A. 

FIG. 4. Model C I-V's for a DC-driven array, as a function 

of D, for an 8 X 8 array. The solid line was obtained with 
D = 0.20, the dashed one with D = 0.99999. /3c = 0.50. 

FIG. 5. FASTCAP model I-V's for a DC-driven array, as a 
function of the ratio of the diameter D to the lattice spacing 
for an 8 X 8 array studied by Whan et a/. [19]. The solid line 
was obtained with D = 0.050, the dashed one with D = 0.970. 
Parameter values axe I3c = 0.50. We use a different diameter 
value than in the two-disk case since they are not exactly the 
same. 

FIG. 6. Model B I-V curves for 40 x 40 array illustrat- 
ing the GCFS of decreasing size and visible at 
t|j, T|j(just). The inset shows the same steps, with the 
t/-coordinate multiplied by 39 to show the precise values of 
the steps. /3c = 0.50, = 0.10, Ao = 1.75. The Ao value 
chosen corresponds to one chosen in previous figures. 

FIG. 7. Four I-V's for Model B, corresponding to 
Ao = 32, 8, 3, 1.75, from left to right, and each displaced for 
clarity successively by 0.5 units to the right. They show how 
the integer n = 1 step is rapidly destroyed as Ao, is reduced 
below the array size (L^ = 16, Lj^ = 16). Just as the integer 
step disappears, a fractional step emerges (a central subject 
of this paper). Parameter values are /3c = 0.5, v = 0.1, and 
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FIG. 8. The n — 1 step width is shown as a function of 
Ao, for Model B. It is clear that once the screening length 
becomes less than half the array size {Lx — Ly — 16), the 
step is destroyed. Values for /3c, i^, and i^c are the same as in 
Fig. |. 

FIG. 9. The width of the fractional |-step is shown as a 
function of u, in Model B. Values for /3c, iac are the same as in 
Fig. and Ao — 1.75. Not shown are values of € [0.3, 0.6], 
for which we have verified that there is no step. 

FIG. 10. The widths of the y-, |-, and |-steps (symbols 
▲ , ■ and ♦ respectively) as a function of Pc, in Model B. 
Parameter values for v, and Ao are the same as in Fig. We 
have also looked for (but not found) other steps but did not 
find them for (5^ < 10.0. 



FIG. 16. We show, as a function of Lx in Model B, the 
value of the screening length which maximizes the width of 
the fractional steps. /3c, iac, and v are as in Fig. |^ The line 
is simply a guide to the eye. 

FIG. 17. The maximal width Ai™c'"' of the fractional steps 
in Model B, is shown as a function of lattice size. The symbols 
• , ■, ♦ represent the first through third steps respectively 
(i.e., 77"^, 77"^) 7~rT etc.). (Not all steps can be obtained 
for all lattice sizes.) This step-width appears to saturate, 
though not at the same lattice size for all steps. 

FIG. 18. Step width A function of disk diameter, 

D, in Model C for a 16 x 16 array. Data are for ^-step 
(the lowest-order step seen). Parameter values are /3c = 0.50, 
V = 0.10. 



FIG. 11. Model B widths of the i- and |-steps (symbols 
• and ▲ respectively) shown as functions of Aq. Values for 
/3c and v are the same as in Fig. ^ We also have looked for 
other steps for all values of Ao < 32. The inset shows an 
I-V showing the hysteresis of the fractional steps in Model 
B. The dashed (solid) line indicates data taken as the cur- 
rent is increased (decreased), as shown by the arrows. Here 
/3c = 0.50, f = 0.10, Ao = 1.75, Suc = 0.005. 



FIG. 12. The width of the f-step as a function of the di- 
ameter, D, in Model C. See Fig. ^ for an explanation of the 
symbol D. Pc = 0.5, v — 0.1. Note that the step exists only 
for a small range of _D-values. 



FIG. 13. The width of the |-step (•) as a function of Ao, 
for Model C with D — 0.994 (this value gives the maximum 
step- width for Ao = 5.66). The step- width in Model B is 
shown for comparison (♦). /3c — 0.50, u — 0.10, Lx ~ Ly = 8. 



FIG. 14. Model C I-V's for individual rows of array show- 
ing |-step. The solid line represents the array average, 
the dashed line represents the edge junctions, while the 
dashed-dotted line represents interior junctions. Parameter 
values are /3c — 0.50, v = 0.10. The inset shows the power 
spectrum for the ^-step in Model B. Parameter values are 
/3c — 0.50, v — 0.10. The drive frequency and harmonics are 
clearly visible, as are other peaks at roughly ^ and We 
have verified that these latter are transients. 



FIG. 15. This figure shows the minimum value of the 
screening length Ac, as a function of Lx, for which the first 
step is visible, in Model B. The data points are shown as •, 
fitted to the dashed line by logAo ~ log LS '*^**'' ''^^ Values 
for /3c, and u axe as shown in Fig. M. 
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